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Using computer simulations of an atomistic glass-forming liquid, we investigate the fluctuations 
of the overlap between a fluid conhguration and a quenched reference system. We find that large 
fluctuations of the overlap develop as temperature decreases, consistent with the existence of the ran¬ 
dom critical point that is predicted by effective held theories. We discuss the scaling of huctuations 
near the presumed critical point, comparing the observed behaviour with that of the random-held 
Ising model. We argue that this critical point directly reveals the existence of an interfacial tension 
between amorphous metastable states, a quantity relevant both for equilibrium relaxation and for 
nonequilibrium melting of stable glass conhgurations. 
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In the search for a theory of the glass transition, an 
important question is whether the slow dynamics of su¬ 
percooled liquids can be explained in terms of the tem¬ 
perature evolution of ‘universal’ thermodynamic quanti¬ 
ties, independent on the details of the material [i|. Two- 
body density correlators, which are central to the theory 
of simple liquids 0, fail for viscous liquids Q. In the 
framework of the random first-order transition (RFOT) 
theory Q, slow dynamics is described in terms of infre¬ 
quent transitions between amorphous free energy min¬ 
ima, separated by large barriers. Structural relaxation 
is interpreted as the nucleation of one metastable state 
into another [^. This physical picture is supported by 
mean-held calculations of a free-energy V{Q). Here, the 
overlap Q measures the similarity between pairs of con¬ 
hgurations: it acts as an order parameter for the glass 
transition Q. This ‘Landau free energy’ describes the 
overlap huctuations and links statics to dynamics. Be¬ 
cause they embody high-order density correlations, over¬ 
lap huctuations are key candidates to construct a ther¬ 
modynamic theory of the glass transition and, a.s such, 
are currently the focus of a large interest |7Hl3l|. The 
central challenge, tackled here, is to understand whether 
the rigorous results obtained in the mean-held limit are 
relevant for realistic, hnite dimensional liquids. 


At mean-held level, V{Q) allows for a compact de¬ 
scription of the liquid-glass phase transition occurring at 
the Kauzmann temperature, T/y, for which Q is the rele¬ 
vant order parameter. By introducing additional external 
helds, such as a coupling e to a reference copy of the sys¬ 
tem, the glass transition at (T = TK,e = 0) transforms 
in a hrst-order transition line ending at a critical point 
at (r„e, > 0) i[n|,[il. Because > T^, this critical 
point might be more easily accessible than Tjy. This sce¬ 
nario occurs in some disordered spin models but has not 
been demonstrated in equilibrium calculations on atom¬ 
istic systems (see BEllil for earlier work). Here, we 
present free-energy calculations that provide direct evi¬ 
dence for such a critical point in a realistic glass-former. 


and hnd a critical behaviour consistent with the univer¬ 
sality class of the random held ising model (RFIM), in 
agreement with theory 12, [l^. Further motivation to 
analyse this critical point is that its absence would di¬ 
rectly establish that a Kauzmann transition does not 
occur at any Tjy > 0. (The reverse is not true: a h- 
nite Tc does not imply a hnite Tr-.) The existence of a 
hrst-order transition line at T < is even more signif¬ 
icant for understanding glassy phenomenology because 
it is associated with phase coexistence between equiva¬ 
lent metastable states. In particular, the corresponding 
interfacial tension, T(T), is an important piece of the 


RFOT theory ^[l3|, but even establishing its existence 
has proven difficult in computational studies [?, 18, 1^ . 
Our approach bypasses these problems and paves the way 
for direct quantitative determinations of T. As a hrst 
step, we analyse cases where stable glass states coexist 
with liquid states, with a large barrier between them that 
we attribute to an interfacial free energy. This hnding is 
relevant to the current effort to understand the melting 
process of ultrastable glasses pH - E^ . 


We consider a well-studied binary mixture of Lennard- 
Jones particles [1^, which exhibits glassy dynamics be¬ 
low T « 1.0. The unit of length is the diameter of the 
larger particles, the unit of energy is the pair-interaction 
scale, the Boltzmann constant is unity, and the total 
number of particles is N. The position of particle i is 
Vi, and C = (ri, r 2 ,..., r^r) denotes a configuration of 
the system. Following [6|, we simulate two copies of the 
system. First, a reference configuration Co is drawn at 
random from the equilibrium state of the model at tem¬ 
perature T'. Then, we perform Monte-Carlo (MC) simu¬ 
lations on a second configuration C, which is biased by a 
field e to lie close to Cq. Specifically, our MC method sam¬ 
ples configurations C according to the Boltzmann weight 
for the Hamiltonian E{C) — eQ{C, Cq), where E is the po¬ 
tential energy and Q{C,C') = ©(« “ ki “ 'f'jD is 

the overlap between configurations C and C'. Here, Q{x) 
is the Heaviside function and we take a = 0.3. Efficiently 
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sampling this Hamiltonian is challenging, and we adopt 
the strategy detailed in I, [13 which combines parallel 
tempering with umbrella sampling and reweighting tech¬ 
niques. Using this approach we were able to study system 
sizes N = 150, 256 for temperatures T > 0.55. 

The temperatures T and T' may be different (Cq need 
not even be thermalized), but the case T = T' has spe¬ 
cial importance. For e = 0, C and Cq are independent. 
As e increases, the configurations become more similar 
and their mutual overlap grows. Our goal is to investi¬ 
gate the presence of singularities in overlap fluctuations 
in the (e, T) plane. To this end, a central role is played by 
the (e-dependent) distribution of the overlap, P{Q). This 
distribution is obtained by first performing a thermal av¬ 
erage over C, and then averaging the results over refer¬ 
ence configurations Cq (we typically use 60 independent 
configurations for each T and N). This double average 
is indicated with simple brackets (•). 

We emphasize that the reference configuration Cq is 
fixed (or ‘quenched’) within each simulation, and sam¬ 
pling is restricted to configuration C. This approach is 
relevant for glassy systems, because we expect Cq to be 
representative of equilibrium metastable states at tem¬ 
perature T'. The configuration C, at temperature T, may 
then either occupy the same metastable state (when Q 
is large) or a different one (when Q is low). If phase 
coexistence occurs at T' = T, then it takes place be¬ 
tween equivalent metastable states. The ‘annealed’ case 
in which both configurations fluctuate is qualitatively dif¬ 
ferent because the high-Q configurations are not repre¬ 
sentative equilibrium states [iliiS]. 

In Fig. [T] we show results for T' = T. The isotherms 
{Q){e,T) shown in Fig. [IJa) recall the relationship be¬ 
tween pressure and density in a liquid-vapor system. On 
lowering the temperature, they become increasingly flat. 
This evolution is more pronounced for the larger sys¬ 
tem. Below the critical point, one would observe, in 
the thermodynamic limit, horizontal tie-lines associated 
with phase coexistence. Fig. mb) shows the suscepti¬ 
bility x(e,T) = N ((Q^) — (Q)^), which exhibits a pro¬ 
nounced maximum at a field e*{T) that shifts to lower 
values when T decreases, suggesting a reduction in the 
thermodynamic cost of localizing the system near a ref¬ 
erence configuration (related to the configurational en¬ 
tropy). Interestingly x(e*,T) grows rapidly as T is low¬ 
ered, and also increases strongly with system size. This 
behavior is expected if the system approaches a critical 
point, but can also be obtained if a correlation length 
larger than the system size develops. The evolution of 
e*{T) is shown in Fig.[TKc) in the (e, T) phase diagram. In 
the liquid-vapor analogy, this corresponds to the critical 
isochore and, for temperatures above the critical point, 
to a Widom line [s^. The observed linear decrease of e* 
with T is consistent with mean-field theories. 

Finally, Fig. m) shows the free energy V{Q) = 
{—T/N) log P{Q) measured at various T, for e = 0. It 



FIG. 1: (a) Dependence of the average overlap (Q) on the field 
e along various isotherms for N = 256 (full) and N = 150 
(dashed). The behaviour is reminiscent of isotherms in a 
liquid-vapor system, (b) Total susceptibility x(e, T) for the 
same parameters showing increasing fluctuations as tempera¬ 
ture is reduced and system size is increased, (c) Temperature- 
dependence of the field e* that maximises y, indicating the 
boundary between low-Q and high-Q regimes, (d) Free en¬ 
ergy V(Q) for the same temperatures as in (a, b) for N = 256, 
with a dashed line indicating linear behavior. 


has the usual convex form at high T with a minimum 
near Q « 0. For lower T and intermediate values of Q, 
its curvature decreases markedly, consistent with the in¬ 
creasing variance of Q. For T < Tc, one expects V{Q) 
to acquire a region of zero curvature, due to Maxwell’s 
construction, similar to the behavior shown at the low¬ 
est studied temperature T = 0.55. Whereas the value of 
V{Q) in the high-overlap regime provides a direct esti¬ 
mate of the configurational entropy 0, the postulated 
surface tension T(T) between states cannot be obtained 
directly from V{Q). 

To gain further insight, we analyse the overlap distri¬ 
bution, P{Q), because it can better reveal critical fluc¬ 
tuations and phase coexistence. We show in Fig. [21(a) 
the temperature evolution of these distributions evalu¬ 
ated along the e*(T) line. The distribution is narrow and 
essentially Gaussian at high temperature, T > 0.7. It be¬ 
comes broader, develops a non-Gaussian flat center, and 
eventually becomes weakly bimodal for T = 0.55, N = 
256. Bimodality is more pronounced for the largest sys¬ 
tem. To investigate this A^-dependence, we calculated the 
Binder cumulant ratio B = {{Q — {Q)Y)‘^ / {{Q — {Q)Y)- 
On a first-order transition line, B increases with N, tend¬ 
ing to unity as N ^ oo whereas above it decreases 
with N towards At Tc, 73 is size-independent with 
a limiting value intermediate between 7 and unity. At 
high temperatures, we find 73 « 7 as expected, and for 
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FIG. 2: (a) Overlap distributions P(Q) measured along the 
e*(T) line for N = 256 (full lines); dotted lines are for N = 
150 and T = 0.6 and 0.55. For the lowest temperature, a 
bimodal structure appears, (b) Behaviour of P{Q) for various 
e, increasing from low-e (blue) to e = e* (black) and large-e 
(purple) for N = 256 and T = 0.55. 
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FIG. 3: (a) Sample-to-sample fluctuations of the distribu¬ 
tions P{Q) for T = 0.55, N = 256 evaluated at the average 
e*. Colors indicate whether the samples are predominately in 
high-, intermediate, or low-Q states. As expected for RFIM 
behaviour, most distributions are unimodal, averaging to a 
total bimodal distribution (black), (b) Temperature depen¬ 
dence of the susceptibilities Xdis and xt measured at e*. The 
solid line is the RFIM prediction, Eq. O- 



T = 0.55 and N = (150,256), we get B = (0.48,0.52), 
consistent with this temperature being close to Tc, al¬ 
though our range of system sizes is rather small. For 
T = 0.55, N = 256 we show PiQ) in Fig. (Kb), which 
evolves from a unimodal distribution at low-Q when 
e < e* to a similarly unimodal distribution at high Q 
when e > e*, with an intermediate bimodal distribution 
at e*. This evolution is qualitatively consistent with the 
crossing of a first-order-like transition line in a finite- 
size system, but systematic finite-size scaling is needed 
to draw conclusions about the thermodynamic limit. 

Random critical points in systems with quenched disor¬ 
der have distinct properties from standard liquid-vapor 
transitions 31|. In particular, the system may behave 


differently for different realizations of the disorder [32l| . 
in our case different reference configurations Cq ■ It is thus 
useful to decompose averages, using {■)co for a thermal 
average over C at fixed Cq, and an overbar (■) for a dis¬ 
order average over Cq. We find that the isotherms {Q)co 
have strong sample-to-sample fluctuations. The overlap 
grows more or less abruptly with e for different samples, 
so that the critical field is a distributed quantity. In 
Fig. [3Ka) we show that at low temperature, T = 0.55, 

the distribution of e% is so broad that the distributions 
Co _ 

P{Q) for most samples at the average value e* = 

are dominated by either their small-Q or large-Q peaks. 

Most of these distributions do become bimodal when the 

field is a justed to foi' each sample separately. These 

results strongly suggest that the quenched disorder plays 

a significant role in the overlap fluctuations. 

In mapping the supercooled liquid to the RFIM, the 
configuration C corresponds to an Ising model in a ran¬ 
dom magnetic field, and the specific realisation of the 
disorder corresponds to Cq. In the RFIM, each realisa¬ 
tion of the random field slightly favors either positive or 
negative magnetisation. As a result, each s amp le typi¬ 
cally has only one thermally-populated state 32| , in con- 


strast to the pure case for which both states are equally 
likely. In large systems, this disorder effect dominates 
the critical fluctuations, which changes the universality 
class from Ising to RFIM. 

For quantitative analysis, we decompose fluctuations 

into thermal, xt = -^((Q^)co “ (*3)co)’ disorder 

_ _ 2 

parts, Xdis = N{{Q)1^ - {Q)ca )■ Clearly, x = Xdis + XT- 
Numerical results, shown in Fig. [Sjb), indicate that 
for the largest system the disorder fluctuations increase 
much more rapidly than thermal ones, 

Xdis oc Xt> (1) 


implying that the disorder should eventually completely 
dominate the behaviour of the total susceptibility, x ~ 
Xdis » XT. Physically, Eq. m can be understood by 
assuming that fluctuations are completely dominated by 
the sample fluctuations of the field In that case, 
standard results for ensemble-dependence of fluctuations 
"ssl-SSl yield Xdis ~ Var(e*)xT/T^, where Var(e*) is 


the variance of among samples. The resulting scal¬ 
ing relation o is characteristic of the RFIM universality 
class. It is exact in mean-field analysis [l^, and the de¬ 
viations predicted from nonperturbative renormalization 
treatments in d = 3 are too small to be numerically 
observable js^]. Interestingly, the ratio Xdis/XT provides 
a measure of the (effective) variance of the field, which is 
a dimensionless measure of the strength of the disorder in 
the effective RFIM description of the critical point [l^ . 
It would be interesting to perform similar measurements 
for other models to perform a quantitative comparison of 
effective theories for different liquids [s^. 

The distribution P{Q) allows an estimation of the 
surface tension T(r), provided the first-order transition 
regime characterized by well-separated peaks can be ac¬ 
cessed. Unfortunately, the computational effort required 
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FIG. 4: (a) Comparison of distributions P{Q) and T = 0.55, 
for T' = T and T' = 0.45, at their corresponding fields e*. (b) 
Comparison a.t T = 0.6 between equilibrated Co at T' — 0.6 
and (biased) inactive Cq. In both cases, the more stable states 
are associated with strongly bimodal P{Q), indicating that 
coexistence between these low-energy states and the equilib¬ 
rium fluid incurs a significant interfacial free energy cost. 


to sample P{Q) for T <IC 0.55 is prohibitive. To cir¬ 
cumvent this difficulty, we exploit the flexibility offered 
by the independent sampling of configurations C and Cq. 
When T' ^ T, mean-field theory predicts that the three- 
parameter space (e, T, T') contains a line of RFIM critical 
points. The phase diagram in Fig. [TJc) corresponds to 
the T' = T plane in that space, but fixing for instance 
T' = 0.45 should yield another critical point {ec,Tc). If 
Tc > Tc, then sampling temperatures below Tc should 
be easier. This is demonstrated in Fig. HKa) where P{Q) 
for (T',T) = (0.45,0.55) shows dramatically enhanced 
bimodality as compared to T' = T = 0.55. The natural 
inference is that the reduction in T' has resulted in an in¬ 
crease of the critical temperature, so that T = 0.55 is now 
effectively deeper into the coexistence region. Thus, we 
interpret the free energy minimum in P{Q) as an inter¬ 
facial cost, now corresponding to the spatial coexistence 
of equilibrium states at temperatures T' and T. 


We show in Fig. HKb) a similar comparison for T = 0.6 
where configurations Cq at taken either from equilib¬ 
rium at T' = 0.6, or sampled from the non-equilibrium 
s-ensemble at that same temperature with a biasing 
field towards atypically low dynamical activity [ 23 ,[ 2 l|. 
In terms of inherent structure energies, these inactive 
states represent glasses with very low fictive tempera¬ 
tures, T' ss 0.4 [2l|. We observe again strong bimodality, 
in contrast to the corresponding equilibrium behaviour, 
which has nearly Gaussian fluctuations. The free-energy 
minimum in P{Q) is the interfacial cost between low-Q 
states (typical of the equilibrium fluid at T = 0.6), and 
high-Q states (stable glassy states). 

The results in Fig. |4Ka,b) are significant because they 
directly evidence the existence of a free-energy cost for 
spatial coexistence between amorphous states. For (T' = 
0.45, T = 0.55) we obtained data for N = 150 and 256 
and find that this cost increases with N. Additional sys¬ 


tem sizes are needed to infer a numerical value for T, via 
the scaling AF ^ TL®. The exponent 9 is not known, 
but should obey 6 < d — I, because interfaces can use 
quenched disorder to optimize their geometry and reduce 
their cost 3lj. Physically the existence of an interfacial 
cost for the coexistence of stable glassy configurations 
with the equilibrium fluid, as in Fig. [3Kb), suggests that 
nucleation and growth is the appropriate mechanism to 
interpret the kinetic stability and melting dynamics of 
ultrastable glasses [H, Hi • Systematic studies along the 
present lines should help understanding kinetic stability 
from first-principles. 

We have presented the results of extensive free-energy 
calculations in supercooled liquids to assess the exis¬ 
tence of thermodynamic singularities in constrained su¬ 
percooled liquids, and to explore their consequences. 
Overall, our results are consistent with effective field 
theories predicting the existence of a random critical 
point, associated (below Tc) with phase coexistence be¬ 
tween metastable states. These results are also consis¬ 
tent with static correlations on length scales comparable 


with our largest system size [39|, |40| . Whatever the sta¬ 
tus of the phase transition in the thermodynamic limit, 
we find that significant static fluctuations are present 
at a temperature T < 0.55, where dynamics is glassy. 
Note that if the critical point in this system is indeed 
at Tc ss 0.55, a significant surface tension T only ex¬ 
ists when T 0.55. Thus, the emergence of activated 
dynamics between metastable states might well fall out 
of the dynamic range currently accessible to simulations. 
Our approach generically allows the determination of T 
in any material (including hard spheres), and extends to 
stable glassy states. This offers the potential for quan¬ 
titative analysis of glass stability based on surface and 
bulk free energies. Overall, while our results do not 
speak directly to the mechanism of structural relaxation 
in glasses and supercooled liquids, they provide direct 
evidence that an interfacial free-energy barrier between 
metastable states is relevant both in the ‘melting’ of sta¬ 
ble glasses and in glassy dynamics at equilibrum. 
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